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ABSTRACT 


Optimal  control  theory  suggests  maintaining  an  orbital  altitude  band  for  Low- 
Earth-Orbiting  (LEO)  satellites  using  periodic  thrusting  than  forced  Keplerian  motion,  i  e. 
a  trajectory  obtained  by  thrust-drag  cancellation  Designing  guidance  algorithms  for  orbit 
maintenance  is  complicated  by  the  nonlinearities  associated  with  orbital  motion.  An 
algorithm  developed  previously  using  thrusters  firing  significantly  off  the  direction  of 
motion  successfully  maintains  an  orbital  band,  but  is  very  inefficient.  This  thesis  develops 
two  different  control  strategies  based  on  the  osculating  orbital  parameters,  taking  a 
conservative  approach  to  keeping  within  altitude  limitations.  Thrust  is  in  the  local 
horizontal  plane,  along  the  direction  of  flight.  Single  and  dual  bum  maneuvers  are 
considered  for  various  bandwidths  and  thmster  sizes.  The  dual  bum  strategy  is  somewhat 
close  to  a  Hohmann  transfer.  The  specified  orbital  band  is  generally  maintained,  with 
some  cases  slightly  exceeding  the  upper  limit.  Propellant  consumption  for  both  maneuvers 
is  significantly  better  than  previous  methods.  This  thesis  shows  that  forward  firing 
thmsters  can  be  used  with  osculating  orbital  parameters  to  obtain  efficiencies  within 
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I.  INTRODUCTION 


The  rising  cost  of  placing  a  satellite  in  orbit  is  creating  a  higher  emphasis  on 
minimization  of  non-payload  mass  to  facilitate  more  payload  items,  and  thus  capability,  on 
a  single  vehicle  Reducing  the  mass  of  propellant  required  to  maintain  a  desired  orbit  has 
the  additional  benefits  for  the  program  in  that  it  also  reduces  the  size  and  mass  of  the  tanks 
needed  to  store  it.  Spacecraft  configuration,  orbit  requirements  and  thruster  control  logic 
dictate  propellant  consumption.  The  propellant  mass  required  for  orbital  maintenance  is 
significantly  higher  for  a  Low-Earth  Orbiting  (LEO)  satellite  than  a  geosynchronous  one 
because  of  the  significant  effects  of  atmospheric  drag  on  the  orbit. 

Consider  the  problem  of  maintaining  a  spacecraft  within  a  prescribed  orbital  band 
in  LEO  (see  Figure  1).  The  conceptually  simplest  thruster  control  logic  is  that  of  a  Forced 
Keplerian  Trajectory  (FKT),  where  the  thrust  cancels  drag.  While  this  method  is  easy  to 
visualize  it  is  difficult  to  implement  because  of  technological  restrictions  on  drag 
estimation,  thrust  vectoring  and  thrust  magnitude  adjustments.  Although  Ross  and 
Melton,  using  optimal  control  theory,  show  that  an  FKT  is  not  an  optimal  maneuver  for 
propellant  consumption  [Ref  1],  it  is  a  good  baseline  to  measure  other  methods. 


Figure  1  Orbit  Radial  Bandwidth 

The  Hohmann  transfer  between  the  lower  and  upper  radial  bands  is  another 
method  of  thruster  control.  Two  thruster  burn  sequences  are  used  in  this  scheme  The 
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change  in  velocity  imparted  at  the  lower  altitude  in  the  first  sequence  places  the  satellite  in 
an  elliptic  transfer  orbit  that  has  an  apogee  of  the  maximum  radius,  while  the  one  at  the 
upper  altitude  circularizes  the  orbit.  Both  thrust  sequences  are  in  the  direction  of  flight. 

Gottlieb's  Eccentricity  Intercept  Targeting  and  Guidance  (EITAG)  [Ref  2] 
program  computes  thruster  bum  duration  and  timing  based  on  the  relationship  between 
eccentricity,  e,  and  semi-major  axis,  a,  as  shown  in  Figure  2.  It  integrates  forward  from 
the  current  s-a  pair  and  backward  from  the  desired  c-ci  pair,  A  two  bum  routine  reboosts 
the  satellite  to  the  desired  circular  orbit.  The  first  bum  achieves  the  intersection  point  of 
the  two  curves  The  vehicle  coasts  until  it  reaches  the  required  position  in  this 
intermediate  orbit  that  enables  the  second  bum  to  move  the  satellite  along  the  final  curve, 


Figure  2.  Eccentricity  Versus  Semi-major  Axis 
In  order  to  maintain  a  radial  band  Pauls  [Ref  3]  and  Wilsey  [Ref  4]  developed  a 
single  bum  control  logic  based  on  orbital  radius,  specific  energy  and  a  thmster  cant  angle 
(the  angle  the  thmst  vector  makes  with  respect  to  the  local  horizontal)  While  the  method 
successfully  controls  the  radial  band  it  causes  the  satellite  to  consume  at  least  three  times 
the  propellant  of  the  FKT  benchmark. 

Since  orbital  parameters  are  derived  from  the  radius  and  its  changes,  in  this  thesis  a 
control  strategy  is  developed  using  the  osculating  (instantaneous)  perigee  and  apogee  radii 
for  a  single  bum  with  the  thruster  firing  in  the  local  horizontal  plane  along  the  direction  of 
motion  A  dual  burn  logic,  using  the  single  bum  strategy  as  the  first  bum  and  the  flight 
path  angle  as  the  control  variable  for  the  second  bum,  is  also  explored 


2 


n.  FORMULA  DEVELOPMENT 


A.  DEVELOPMENT  OF  THE  EQUATIONS  OF  MOTION 

Assuming  the  satellite  is  a  non-lifting  (blunt)  body  simplifies  analysis  in  that  drag  is 
the  only  uncontrollable  nonconservative  force  being  considered.  Specifying  planar 
motion,  the  initial  orbit  as  circular,  defining  the  problem  as  a  two-body  problem  and 
neglecting  all  other  possible  perturbations  facilitates  ease  of  formulation.  The  orbital 
dynamics  can  be  considered  in  a  polar  coordinate  system  as  shown  in  Figure  3. 


Figure  3 .  Orbital  Coordinate  System 
The  equations  of  motion  are 


=  (1) 

a,  =  S§  (2) 

where  and  a,  are  the  radial  and  transverse  components  of  the  inertial  acceleration, 
and  ZF^are  the  sums  of  the  respective  forces,  and  m  is  the  satellite  mass.  Drag  and  thrust 
can  be  broken  into  components 


3 


Dr  =  -Dsmy 


(3a) 


Dt  =  -D  cosy  (3b) 

Tr^Tsma  (4a) 

Tt  =  Tcosa  (4b) 

where  y  is  the  flight  path  angle  (the  angle  between  the  transverse  axis  and  the  velocity 
vector)  and  a  is  the  thruster  elevation  angle  (the  angle  between  the  transverse  axis  and  the 
thrust  vector).  Expressing  the  equations  of  motion  in  polar  coordinates  and  substituting  in 
Equations  3  and  4  yields 


a2  _  g  I  Dr  , 

^  ~  ~ ^2  W  W 

Qr  +  2Qr=^  +  i  (6) 


where  g  is  the  geocentric  gravitational  constant  and  0  is  the  angular  position  of  the 
satellite  as  depicted  in  Figure  3, 

B.  NONDIMENSIONALIZATION  OF  THE  EQUATIONS  OF  MOTION 

It  is  useful  to  nondimensionalize  the  equations  of  motion  listed  above  in  order  to 
limit  the  effects  of  computational  errors  when  dealing  with  large  numbers.  This  also 
allows  better  analysis  of  the  effects  of  the  variation  of  parameters. 

1.  Definitions 

Base  units  in  mass,  length  and  time  are  required  to  nondimensionalize  the  problem 
The  base  mass,  /w^,  is  the  initial  mass  of  the  satellite.  The  base  length,  is  the  initial 
semi-major  axis  of  the  actual  orbit.  Since  the  satellite  travels  above  and  below  the  center 
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of  the  desired  radial  band  the  initial  satellite  position  is  at  the  upper  limit  of  the  bandwidth 
The  base  time,  is  the  period  of  the  initial  orbit 

tb  =  271 

(7) 

The  nondimensionalized  mass,  radius  and  time  are  thus 

—  m 

fn  — 

(8) 

(9) 

'‘=i 

(10) 

2.  Equation  Nondimensionalization 

The  satellite's  position  variables  and  their  derivatives  are  nondimensionalized  using 
the  base  units 


^  _  H-p/ 

~  dt  ~  d{tbt)  ~  tbcCt  h 


(H) 


^£L-±(^  -Af Hr/)  -  H-^(r/\  -  =  Hy/l 

~  dt^  dt\dtj  dt\tb  J  ''f  d(tbty  ^  tl  dt  ^ 


(12) 


e  =  e 


(13) 


A d^  _  d^  _  dS _ ^0^ 

~  dt  ~  d{tb't)  ~  tbdi  ~  h 


(14) 


dt'^  ~  dt\dtj 


7 

dt\tb 


1  d 
hd{tbt) 


t\  eft 


-f 


(15) 
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where  primes  represent  a  differentiation  with  respect  to  nondimensionalized  time 
Substituting  these  equations  and  Equations  3  and  4  into  the  equations  of  motion  yields 


3.  Nondimensionalization  of  Nonconservative  Forces 

The  nonconservative  forces  (thrust  and  drag)  are  to  be  nondimensionalized  by 
multiplying  and  dividing  by  the  base  units  required  to  remove  the  dimensions 
a.  Nondimensionalized  Thrust 


Nondimensionalized  thrust  is  defined  by 


T-f  _  *  0 


b.  Nondimensionalized  Drag 
Nondimensionalized  drag  is  given  as 


D  =  t^D 


4.  Nondimensionalization  of  Other  Parameters 

The  drag  the  satellite  encounters  depends  solely  on  the  atmosphere  the  orbit  passes 
through.  Drag  is  defined  by 

D  =  \pACDy^  (24) 

where  p  is  the  atmospheric  density,  A  is  the  effective  surface  area  of  the  vehicle,  Q  is  the 
drag  coefficient  associated  with  the  spacecraft,  and  v  is  the  instantaneous  velocity.  An 
exponentially  decaying  density  is  a  common  atmospheric  model  and  is  given  by 

p  =  (25) 

where  is  the  density  at  the  reference  altitude,  and  H\s  the  atmospheric  scale  height. 
Nondimensionalized  scale  height  is  defined  as 


A  base  density  is  defined  as 


9h  = 


7 


The  nondimensionalized  density  is 


PlQ-^k^-^ref) 


The  ballistic  coefficient  is  defined  as 


B  = 


mb 

ACd 


thus  the  nondimensionalized  ballistic  coefficient  is 

^  “  9bn  ACdpbrb 


The  spacecraft's  velocity  is 


=r^  +r^  0^ 


and  is  nondimensionalized  using  Equations  9,  1 1  and  14  to 

The  nondimensionalized  thrust-to-drag  ratio  is 

J_ 

Do 

where  Do  is  the  nondimensionalized  drag  for  the  FKT  orbit  The  change 
by 


(28) 

(29) 

(30) 

(31) 

.2 

jv^  (32) 

b 

(33) 

in  mass  is  given 
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W=  -1 


Nondimensionalizing  yields 


—  /  tb 

m' 


T  Tmtrb  tb _ Trb 

hpg  Tt] 


5.  Nondimensionalized  Equations 

Substituting  Equations  22  through  32  into  the  equations  of  motion  produces 


r"  =  g'^P- 


^  iy^siny  fsina 

r  )  2mB  ^ 


r 


^  ^v^cosy  Teas  a 

ImrB 
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m.  COMPUTER  MODEL  DEVELOPMENT 


A.  ORIGINAL  development 

The  original  program  was  developed  by  Pauls  [Ref.  3],  and  modified  by  Wilsey 
[Ref  4]  and  Gardner  [Ref  5],  It  was  FORTRAN  based,  and  simulated  the  motion  of  a 
spacecraft  in  orbit  by  utilizing  a  fourth-order  Runge-Kutta  numerical  integration  routine. 
Wilsey  [Ref  4]  nondimensionalized  the  equations  and  used  four  state  variables  to  reduce 
the  two  second  order  equations  of  motion  to  four  first  order  equations;  Gardner  [Ref.  5] 
redefined  the  variables  and  added  the  mass  to  the  state  variables  and  the  mass  equation  as 
a  first  order  differential  integrated  by  the  Runge-Kutta  routine. 

B.  STATE  VARIABLE  DEFINITIONS 

Maintaining  a  four  element  state  vector  based  on  the  position  variables  and  leaving 
the  mass  as  a  separate  issue,  the  state  vector  is  defined  as 


XI  =r  (38) 

X2  =r^X\  (39) 

X3  =  e  (40) 

X4  =e=JC3  (41) 

Their  nondimensionalized  counterparts  are 

Xi=r  (42) 

X2=r'  =  x\  (43) 

X3  =  0  (44) 
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JC4  =  0^  =  ^3 


C.  NONDIMENSIONALIZED  STATE  VARIABLE  EQUATIONS  OF 
MOTION 

Rewriting  the  equations  of  motion  (Equations  38  and  39)  by  substituting  in  the 
nondimensionalized  state  variables.  Equations  42  through  45,  yields 


X2  =X\xl 


^^v^siny  ^  Tsina 
2mB  ^ 


-<  = 


2x4X2 


Icosa 

2mxiB 


The  flight  path  angle,  y,  is  dependent  on  the  radial  and  tangential  components  of  the 
velocity  vector  and  can  be  expressed  in  terms  of  the  state  variables 


X2 

smy  =  - 


X1X4 

cosY  =  ^ 


Thus  the  equations  appear  in  the  computer  program  as 


Xi  =  X2 


-/  —  —2  47c^ 

X2  ~  2 


2mB 
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(52) 


_/  2X4X2 

^4  -  3ci 


2mB 


+ 


T  cosa 
lm\ 


(53) 


D.  PROGRAM  ORGANIZATION 

The  program,  written  in  MATLAB,  is  broken  into  several  functional  subroutines  to 
provide  a  fluid  and  logical  flow  of  data  through  required  calculations.  Initialization  of 
constants  and  variables  is  followed  by  the  computation  of  drag  at  the  Keplerian  orbit 
radius  and  then  the  main  program.  Nondimensionalized  velocity  is  computed  in  LEOVEL 
using  the  radial  and  tangential  velocity  components.  The  exponential  term  of  drag  is 
evaluated  in  LEODRAG.  LEOLOCAL  determines  the  individual  pieces  of  the  equations 
of  motion  for  easy  assimilation  in  LEOXDOT,  where  the  derivatives  of  the 
nondimensionalized  state  variables  are  calculated.  The  Runge-Kutta  is  completed  in 
LEORK4A,  LEORK4B,  LEORK4C,  and  LEORK4D,  with  the  program  recycling  through 
from  LEOVEL.  After  the  RK4  routines  have  produced  the  new  nondimensionalized  state 
variables  the  propellant  consumption  for  that  step  is  computed,  then  the  osculating  orbital 
parameters  are  determined  by  LEOPARMl  and  LEOPARM2.  The  thruster  firing  logic 
routine  in  effect  is  called  (either  LEOTFLPH  alone  or  in  combination  Avith  LEOTFLAP 
for  the  dual  bum  scenario)  to  determine  whether  or  not  the  thmster  is  active  during  the 
upcoming  step.  A  flowchart  for  the  program  is  shown  in  Figure  4. 

E.  PROGRAM  VALIDATION 

All  aspects  of  the  code  were  verified  since  the  program  was  written  in  a  different 
language  from  the  previous  versions.  The  process  Wilsey  [Ref  4]  utilized  was  duplicated. 

1.  Initial  Validation  with  No  External  Forces 

With  no  atmospheric  affects  to  hinder  orbital  motion  the  specific  energy  and 
angular  momentum  for  both  elliptic  and  circular  orbits  must  remain  constant.  Results  for 
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these  values  are  shown  in  Figures  5  and  6  for  the  elliptic  orbit,  and  Figures  7  and  8  for  the 
circular  orbit  Small  fluctuations  are  impossible  to  see  on  the  plots,  so  closer  analysis  is 
called  for  Comparison  of  the  output  of  the  program  to  the  initial  values  for  these 
arguments,  depicted  in  Figures  9  and  1 0  for  the  elliptic  case  and  Figures  1 1  and  1 2  for  the 


circular  one,  shows  minor  fluctuations  on  the  order  of  lO"'^  with  respect  to  the  base  for  the 
elliptic  orbit  and  zero  for  the  circular  orbit.  The  fluctuations  follow  no  discernible  pattern 
and  are  of  a  relative  magnitude  that  is  known  to  be  beneath  MATLAB’s  computational 
limits  and  are  viewed  by  the  processor  as  near  zero.  This  suggests  they  are  computational 
noise  generated  by  the  numerous  calculations  of  infinitesimal  numbers,  and  the  main  body 
of  the  program  is  validated. 
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Figure  7  Specific  Energy  of  a  Drag  Free  Circular  Orbit 
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2.  Validation  of  Drag 

The  program  was  run  using  a  constant  atmospheric  density  model  The  change  in 
semi-major  axis  and  velocity  were  compared  with  analytic  values  obtained  by  using 
equations  from  Wertz  [Ref.  6]  that  approximate  the  changes  for  a  satellite  in  a  circular 
orbit  experiencing  drag 


Aa  =  -2%^Qa^  (54) 

Av  =  K^pav  (55) 

The  equations  were  modified  to  facilitate  ease  of  computation  by  making  use  of  the 
definitions  of  the  ballistic  coefficient  and  its  nondimensionalized  counterpart,  Equations  29 
and  30,  and  noting  that  the  mass  and  density  are  constant  and  that  the  initial  semi-major 
axis  is  the  base  radius,  to  produce 


Aa  = -27t-srp»a  = -2"| 


(56) 


.  ACd  Wj,  B  _v 

Av  =  ;c— = 


(57) 


A  ten  orbit  test  produced  the  computational  values,  the  analytic  values  were  calculated 
based  on  the  computational  semi-major  axis  and  velocity  at  the  start  of  each  orbit.  Percent 
differences  were  computed  for  both  the  change  and  the  element  The  results  of  the  test, 
shown  in  T able  1 ,  show  that  the  percent  difference  in  all  cases  is  less  than  six-tenths  of  one 
percent,  and  the  drag  routine  is  validated. 


19 


ORBIT 

%  ASMA 
DIFFERENCE 

%SMA 

DIFFERENCE 

%AV 

DIFFERENCE 

%V 

DIFFERENCE 

1 

-0.0063 

-0.0158 

0.0251 

0.0316 

2 

0.0063 

0.0158 

0.0881 

0.1105 

3 

0.0188 

0.0474 

0  1511 

0.1895 

4 

0,0314 

0.079 

0.2141 

0.2685 

5 

0.0439 

0  1106 

0.2772 

0.3474 

6 

00565 

0,1422 

0.3404 

0.4261 

7 

0,0696 

0.1739 

0.4033 

0,5047 

8 

0.0817 

0.2055 

0.46^ 

0.5829 

9 

0.0942 

0.2371 

0.5285 

0.6605 

10 

0.1068 

0.2688 

0,5903 

0.7373 

Table  1  Percent  Difference  in  Change  and  Values  for  SMA  and  Velocity 
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rV.  DEVELOPMENT  OF  THE  ORBIT  MAINTENANCE  PROCEDURE 


Ross  and  Melton  [Ref.  1]  demonstrated  the  non-optimality  of  a  Forced  Keplerian 
Trajectory  with  respect  to  propellant  consumption.  Pauls  [Ref  3]  and  Wilsey  [Ref  4] 
showed  the  adequacy,  but  non-optimality  of  a  single,  off  axis  bum  controlled  by  radius  and 
specific  energy  A  different  method  of  control  must  be  considered  in  order  to  approach 
optimality  while  maintaining  the  required  orbital  band  A  review  of  orbital  mechanics  and 
Wilsey’s  [Ref  4]  plots  for  radius  and  eccentricity  suggested  a  potential  solution. 

A.  THRUSTER  CANT  ANGLE 

A  firing  thruster  imparts  energy  to  the  satellite  at  the  expense  of  propellant,  thus 
maximizing  the  energy  gained  while  minimizing  the  bum  time  will  reduce  propellant 
consumption.  The  energy  equation  for  an  orbital  body  is 


\l 

8  —  ^  r 


The  second  term  of  the  equation  dominates,  and  orbital  energy  is  negative.  As  the  radius 
decreases  the  second  term  grows,  and  the  energy  increases  negatively.  Increasing  the 
radius  decreases  the  magnitude  of  the  energy.  Zeleny  [Ref  7]  notes  that  Equation  58  can 
be  expressed  as 


8  =  ^  V  •  V  -7 


The  change  in  specific  energy  with  respect  to  time  is 


—  =  V  • 

dt 
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where  ^  is  the  angle  between  the  acceleration  vector  and  the  tangential  axis.  Thrust  is 
larger  than  drag  because  a  thrust  equal  to  drag  would  require  a  FKT  to  maintain  orbit. 
Figure  13  shows  the  results  at  two  different  thruster  cant  angles.  If  the  cant  angle  is  small, 
the  velocity  will  not  move  radically  above  the  tangential  axis,  and  the  second  term  of 


Figure  13.  Resulting  Velocities  at  Two  Thruster  Cant  Angles 


Equation  60  will  be  very  small  in  comparison  with  the  first  term.  If  the  cant  angle  is  high 
the  initial  response  will  be  reduced,  since  both  terms  are  small  A  sufficiently  long  bum 
would  place  the  velocity  and  the  thrust  almost  in  the  same  direction,  but  the  objective  is  to 
reduce  bum  times.  The  maximum  change  of  orbital  energy  is  obtained  when  the  velocity 
and  acceleration  are  nearly  parallel.  For  decay  of  a  circular  orbit  due  to  drag  the  flight 
path  angle  remains  small,  between  0"  and  -T,  provided  the  initial  altitude,  bandwidth 
(separation  between  maximum  and  miniumum  radii)  and  vehicle  configuration  are 
sufficient  to  preclude  extreme  affects.  If  the  path  becomes  elliptic,  or  the  thruster  is  small 
and  the  bandwidth  larger  than  the  thruster  can  achieve  in  one  orbit,  the  flight  path  angle 
may  vary  between  T  and  -1“  Setting  the  thmster  cant  angle,  a,  equal  to  zero  enables  the 
satellite  to  sufficiently  handle  both  positive  and  negative  flight  path  angles  wdth  equal 
ability,  more  so  than  setting  it  1®  to  either  side.  To  achieve  the  same  change  in  energy  a 
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thruster  firing  at  a  significantly  larger  angle  would  have  to  fire  longer,  consuming  more 
propellant.  The  thruster  cant  angle  is  therefore  set  at  zero. 

B.  ECCENTRICITY  CONTROL 

The  imposition  of  a  maximum  and  minimum  radius  limits  the  eccentricity  the  orbit 
may  obtain  without  having  a  portion  of  the  orbit  exceed  either  boundary.  If  the 
eccentricity  surpasses  the  limit,  the  thruster  control  logic,  whatever  it  may  be,  will  be 
forced  to  fire  the  thrusters  more  often  than  if  the  eccentricity  remains  within  orbital 
bandwidth  requirements.  The  limiting  eccentricity  is  given  by 

_  max- min  .. , .. 

~  Rmax+Ram  ^  ’ 

For  a  thruster  control  logic  to  be  remotely  close  to  optimal  it  must  avoid  allowing  the 
osculating  eccentricity  exceeding  the  maximum  value. 

C.  SINGLE  BURN  THRUST  CONTROL  LOGIC 

Since  the  orbital  parameters  of  interest  are  based  on  the  radius,  orbital  angle  and 
their  respective  rates  of  change,  and  since  the  desired  outcome  is  the  maintenance  of  a 
prespecified  orbital  band  by  limiting  the  eccentricity,  the  control  logic  utilizes  these  values 
to  obtain  the  desired  results. 

1.  Thruster  Firing  Criteria 

Acknowledging  that  the  satellite  would  descend  beneath  the  minimum  radius  if  the 
thruster  was  not  fired  prior  to  reaching  it  dictates  a  bum  start  before  the  lower  radial  limit 
is  reached.  Atmospheric  density  profile,  orbital  bandwidth,  thmster  capacity  and  cant 
angle,  and  satellite  configuration  all  play  a  role  in  how  far  the  vehicle  drops  from  thmster 
initiation  until  the  radius  ceases  to  decrease.  Starting  the  bum  when  the  osculating  perigee 
radius  drops  below  the  minimum  radius  will  at  worst  delay  the  time  until  the  satellite 
passes  the  minimum  value  or  cause  the  satellite  to  begin  to  climb  well  before  the  lower 
radial  limit  is  reached.  Any  result  in  between  the  two  is  reasonably  acceptable  and  the 
thmster  firing  criteria  is  established. 
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2.  Thruster  Turn  Off  Criteria 

Gottlieb  [Ref  2]  shows  the  relationship  between  semi-major  axis,  a,  and 
eccentricity,  e,  during  thrusted  flight.  If  the  thruster  is  allowed  to  fire  until  either  portion 
of  the  a-e  pair  is  beyond  its  limiting  value  the  satellite  will  in  all  likelihood  either  exceed 
the  bandwidth  restrictions  or  require  excessive  bums  to  maintain  it.  Turning  the  thruster 
off  when  the  osculating  apogee  radius  exceeds  the  maximum  radius  implies  that  the  orbital 
radius  will  never  exceed  the  maximum  value.  This  is  mainly  because  the  effect  drag  will 
have  on  reducing  the  apogee  radius,  and  partly  due  to  the  change  in  apogee  for  any  given 
instant  being  small  enough  so  that  the  osculating  apogee  does  not  significantly  exceed  the 
limit.  The  error  in  overshooting  the  upper  limit  is  small  considering  the  second  factor,  and 
the  thruster  turn  off  criteria  is  established.  Figure  14  is  the  logic  flowchart  for  the  single 
bum  maneuver. 

3.  Resulting  Orbit 

Since  the  thrusts  occur  based  on  the  osculating  perigee  falling  below  the  minimum 
limit  the  satellite  may  not  be  at  a  perigee  position  when  the  thmster  is  activated. 
Additionally,  the  velocity  change  will  not  necessarily  correspond  to  that  of  a  Hohmann 
transfer.  Changes  in  the  perigee  radius  and  apogee  radius  may  occur  simultaneously. 
Raising  the  instantaneous  radius  by  increasing  the  apogee  radius  is  the  primary  intent  of 
the  bum,  but  changes  in  the  perigee  radius  are  also  important.  Raising  the  perigee  radius 
along  with  the  apogee  radius  places  the  satellite  in  a  more  desirable  orbit  because  it 
increases  the  time  it  takes  for  the  perigee  radius  to  decay  below  the  minimum  limit, 
increasing  the  period  between  bums 

The  change  in  perigee  radius  can  be  determined  from  analysis  of  the  semi-major 
axis,  eccentricity,  angular  momentum  and  specific  energy.  Zeleny  [Ref  7]  states  the 
change  in  eccentricity  and  angular  momentum  from  a  change  in  velocity  are 

de  =  +  Izhdh)  (62) 
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where  h  is  the  angular  momentum.  The  angular  momentum  is 


Figure  14.  Single  Bum  Maneuver  Thruster  Control  Flow 


and  thus 

/z  =  rv  sin  f  I  -  y)  =rvcosy 


The  change  in  angular  momentum  can  be  expressed  as 


dh  =  ra  cos  ^ 


(66) 


Substituting  Equations  58,  60,  65  and  66  into  Equation  62  yields 

2 

de  =  -^[rv^a  cos  y  cos  (y  -  ^  ^  cos  y 

+rv^a  cos  ^  -  l\xa  cos  ^  (67) 

Since  all  the  angles  in  this  equation  are  zero  or  approximately  zero,  and  the  terms  outside 
the  brackets  are  positive,  the  sign  of  the  equation  can  be  determined  from 

sign(i/e)  =sign  2a(rv^  -  (x)  +  ^  (68) 

Zeleny  [Ref  7]  states  that  the  change  in  semi-major  axis  is 

da  =  -^dz  (69) 

28^ 

The  perigee  radius  is  defined  as 

rp  =  a(\  -  e)  (70) 

so  its  change  is 

drp  =  (1  -  e)da  -  ade  (7 1 ) 

From  Equation  68  the  change  in  eccentricity  varies  from  positive  to  negative  As  the 
velocity  increases  with  a  bum,  the  change  in  eccentricity  is  more  likely  to  become  positive. 
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and  attempt  to  reduce  the  change  in  perigee  radius  Eccentricity  is  generally  very  small, 
and  changes  to  it  even  smaller  for  any  given  instant  in  time.  The  change  in  semi-major 
axis  is  not  small,  in  comparison.  The  first  term  in  Equation  71  is  on  the  order  of  the 
change  in  semi-major  axis  while  the  second  term  is  significantly  smaller  The  overall  result 
is  a  positive  change  in  the  perigee  radius,  placing  the  satellite  in  a  much  more  desirable, 
but  still  eccentric,  orbit. 

D.  DUAL  BURN  CONTROL  LOGIC 

The  single  bum  strategy  described  above  leaves  the  satellite  in  an  elliptic  orbit.  If 
the  perigee  radius  has  not  been  sufficiently  raised,  drag  will  soon  lower  it  beneath  the 
minimum  limit  and  the  thruster  will  fire  again.  Reducing  the  eccentricity  by  conducting  a 
second  bum  as  the  satellite  approaches  apogee  could  result  in  an  increase  in  the  period  of 
radius  boosting  bums,  potentially  saving  propellant  The  objective  of  the  second  bum  is  to 
decrease  eccentricity,  but  not  necessarily  make  the  orbit  circular 

1.  Thruster  Firing  Criteria 

As  the  satellite  approaches  apogee  the  thruster  must  fire  to  increase  both  the 
semi-major  axis  and  the  perigee  radius.  The  flight  path  angle,  y,  is  positive  during  the 
perigee  to  apogee  transition,  and  decreases  as  the  vehicle  nears  apogee.  Waiting  for  the 
flight  path  angle  to  become  negative  before  initiating  thrust  could  raise  the  apogee  radius, 
placing  the  satellite  on  a  path  that  exceeds  the  maximum  radial  limit.  Firing  the  thruster 
when  it  drops  below  a  minimum  positive  value  seems  to  be  the  pmdent  course  of  action. 
An  arbitrary  value  was  selected  for  the  firing  criteria,  and  adjusted  until  it  minimized 
eccentricity  for  the  smallest  thruster-bandwidth  combination.  Using  the  resulting  firing 
angle  with  larger  thrusters  and  bands  yielded  higher  than  desired  eccentricities.  More 
powerfijl  thmsters  change  the  path  quicker,  and  should  require  less  time  to  reduce  the 
eccentricity.  Larger  bandwidths  increase  the  maximum  allowable  eccentricity,  producing 
boost  orbits  with  a  correspondingly  lower  velocity  at  apogee.  Longer  bums  are  required 
at  the  smaller  velocity  to  raise  it  to  a  value  close  to  that  of  a  circular  orbit  at  the  apogee 
radius.  Modifying  the  base  flight  path  firing  angle  by  a  ratio  of  the  bandwidth-to-thrust 


27 


ratio  proved  insufficient.  Using  the  square  of  the  bandwidth  in  the  numerator  of  the  ratio 
produced  more  favorable,  but  not  perfect,  results. 

2.  Thruster  Turn  Off  Criteria 

If  the  thruster  fires  too  long  the  apogee  radius  will  rise,  and  the  path  potentially 
exceed  the  upper  radial  limits.  If  the  bum  is  too  short  the  perigee  radius  will  not  be 
increased  sufiBciently  to  affect  its  decay-to-minimum  time  Controlling  the  end  of  the  bum 
by  comparing  the  perigee  radius  to  the  radius  is  insufficient  because  it  does  not  consider 
the  effects  on  the  apogee  radius,  and  the  bum  might  be  long.  Examination  of  the  flight 
path  path  angle  proves  more  interesting. 

The  tangent  of  the  flight  path  angle  is  given  by 

t:,  =  Xmy  =  f-  (72) 

Its  rate  of  change  is 

dr  r  dr  rdQ  1 

=  7i-76-7e  *  * 

Near  apogee,  with  no  thrust,  the  values  of  the  variables  inside  the  expression  are 


r  »  0 

(74a) 

(74b) 

r<  0 

(74c) 

0>O 

(74d) 

0<O 

(74e) 
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and  the  change  in  the  tangent  of  the  flight  path  angle  is  negative  The  thruster,  as  seen 
from  the  equations  of  motion,  changes  the  rates  in  Equations  74c  and  74e  to  positive 
values.  Since  radial  velocity  is  approximately  zero,  the  radial  acceleration  changes  govern 

Equation  73,  eventually  causing  c/^to  become  positive, 

Continuous  thrusting  at  the  chosen  cant  angle  increases  both  r  and  v.  At  apogee 
thrusting  has  more  affect  on  velocity  than  radius  The  result  is  that  the  orbit  tends  to 
circularize,  then,  as  shown  previously,  the  semi-major  axis  increases,  and  the  orbit 
becomes  elliptic  again.  This  ties  in  with  the  change  the  flight  path  angle  undergoes,  and 
thus  when  it  begins  to  increase,  the  eccentricity  has  reached  a  minimum.  The  thruster  cut 
off  criteria  is  determined.  Figure  1 5  depicts  the  flowpath  for  the  second  bum  in  the  dual 
bum  thmster  control  logic. 

3.  Resulting  Orbit 

The  resulting  orbit  will  be  nearly  circular  and  at  or  near  the  maximum  radial  limit 
Since  the  perigee  radius  is  almost  equal  to  the  apogee  radius,  the  decay  pattern  should 
allow  for  more  time  between  boosting  bums. 
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Figure  15.  Dual  Bum  Strategy  Second  Bum  Logic 
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V.  ANALYSIS  AND  RESULTS 


The  satellite  is  initially  positioned  at  the  upper  radius  limit,  as  though  a  launch 
vehicle  had  placed  it  there.  FKT  analysis  is  done  at  the  lower  limit,  midband,  and  the 
initial  radius.  Configuration  changes  are  limited  to  the  thruster  sizing  and  the  limited 
effects  on  nondimensionalization  caused  by  the  different  starting  altitudes.  The  following 
parameters  were  used  in  all  cases. 


Midband  Radius: 

R^  +  260  km 

Cant  Angle: 

0“ 

Ballistic  Coefficient: 

150  kg/m^ 

Scale  Height: 

46.9  km 

Specific  Impulse: 

300  s 

Initial  Mass: 

20,000  kg 

Base  Density  Altitude.  250  km 

Base  Density: 

7.248  X  lO'"  kg/m" 

Density  Factor: 

12.47 

Bandwidths  were  selected  as  2,  10,  25,  and  75  km  (0.3  x  10^,  1.5  x  10  ,  3.8  x  10  ,  and 
11.2  X  10'^  distance  units,  respectively).  Thruster  sizes  were  limited  to  40,  80,  160,  320, 
640,  and  1280  N  Table  2  shows  the  nondimensionalized  counterparts  to  thruster  size. 
The  program  is  run  for  100  orbits  for  the  first  three  bandwidths,  and  200  orbits  for  the  last 
bandwidth.  The  equations  of  motion  are  updated  1 000  times  and  the  output  sampled  ten 
times  per  orbit 

A.  SINGLE  BURN  STRATEGY 
L  Band  Maintenance 

Pauls  [Ref  3]  and  Wilsey  [Ref  4]  demonstrated  that  their  single  burn  control  logic 
did  not  permit  forward-firing  thrusters  to  maintain  an  orbital  band.  The  ability  of  the 
selected  Thruster  Firing  Control  Logic  (TFCL)  to  keep  the  satellite  in  a  bandwidth  at  all 
must  first,  therefore,  be  established.  Initial  runs  demonstrated  that  the  logic  successfully 
maintained  the  radius  between  two  boundaries  Having  shown  the  capacity  to  constrain 
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the  position,  the  TFCL  had  to  be  evaluated  on  its  ability  to  hold  specified  limits. 
Variations  in  bandwidth  and  thrust  magnitude  were  applied.  Figures  16  through  39  show 
the  radius  versus  orbit.  The  TFCL  kept  the  vehicle  at  or  above  the  minimum  radius 
regardless  of  the  band  and  thruster  size  selected.  Using  thrusters  smaller  than  40  N  could 
cause  the  path  to  drop  under  the  limit.  Staying  beneath  the  upper  radial  limit  depended  on 
the  thrust  and  the  bandwidth,  as  well  as  the  time  step  size.  High  thrust  at  small 
bandwidths  caused  the  radius  to  exceed  the  maximum  boundary.  Table  3  summarizes  the 
single  bum  band  maintenance  results. 


Bandwidth 

[km] 

Thmst  [N] 

40 

80 

160 

320 

640 

1,280 

Nondimensionalized  Thmst 

2 

8.73 

17.46 

34.92 

69.85 

139,69 

279.39 

10 

8.74 

17.48 

34.97 

139.87 

279.74 

25 

8.76 

280.37 

75 

8.83 

35.31 

70.62 

141.24 

282.48 

Table  2.  Nondimentionalized  Thrust  by  Bandwidth 


Bandwidth 

[km] 

Thmst  [N] 

40 

80 

160 

320 

640 

1,280 

2 

SAT 

SAT 

SAT 

+  10% 

+10% 

+51% 

10 

SAT 

SAT 

SAT 

SAT 

+4% 

+18% 

25 

SAT 

SAT 

SAT 

SAT 

SAT 

+2% 

75 

SAT 

SAT 

SAT 

SAT 

SAT 

SAT 

SAT  =  Satisfactory  Band  +_%  =  _  %  Above  Maximum  Radius  Limit 

Maintenance  (Percent  of  Bandwith)  _ 


Table  3  Single  Bum  Strategy  Band  Maintenance  Summary 

2.  Orbital  Path 

The  trajectory  caused  by  the  single  bum,  as  expected,  is  generally  elliptic.  The 
orbital  radius  follows  no  immediately  discernible  pattern,  and  neither  does  the  bum 
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pattern.  In  some  instances  TFCL  causes  the  thruster  to  fire  as  the  vehicle  approaches 
apogee,  inadvertently  reducing  the  eccentricity  of  the  orbit,  as  demonstrated  in  Figures  40 
and  41.  This  is  because  the  routine  is  only  looking  at  the  osculating  perigee  radius  for 
firing  criteria,  not  the  vehicle's  position  relative  to  apogee  or  perigee.  Since  the  orbit  is 
not  permitted  to  exceed  the  maximum  eccentricity  imposed  by  the  radial  limits,  bums  are 
required  less  frequently  than  by  the  Pauls-Wilsey  method 

The  method  utilizes  the  full  bandwidth  during  the  course  of  the  run,  but  does  not 
use  all  available  distance  units  on  every  opportunity.  In  the  smaUer  bands,  particularly  the 
2  km  (0.3  X  10'^  distance  units),  the  logic  uses  as  little  as  half  the  available  bandwidth.  In 
cases  with  a  larger  minimum-maximum  separation,  approximately  90%  of  the  band  is  used 
when  bums  occur. 

B.  DUAL  BURN  STRATEGY 

1.  Band  Maintenance 

Again  the  logic  had  to  be  tested  to  ensure  it  was  capable  of  maintaining  a  band. 
Runs  showed  the  dual  bum  maneuver  kept  the  vehicle  within  radial  limits  The  new 
strategy  was  put  through  the  same  bandwidth-thmster  variation  sequence  the  single  bum 
control  method  was.  Plots  of  the  resulting  radius  versus  orbit  are  shown  in  Figures  42 
through  65,  Since  the  single  bum  strategy  was  used  as  the  initial  bum,  the  radius  was  not 
expected  to,  and  did  not,  go  below  the  minimum.  The  timing  of  the  second  bum  did  not 
always  meet  the  needs  of  the  particular  osculating  orbital  elements,  and  periodically  led  to 
the  upper  radial  limit  being  exceeded.  Review  of  the  bum  initiation  time,  the  flight  path 
angle  and  the  satellite's  position  with  respect  to  apogee  revealed  thmsting  began  earlier 
than  required.  Thmst  initiation  for  the  second  bum  is  based  on  a  scaled  flight  path  angle. 
The  scaling  ratio  was  too  large,  indicating  the  method  favors  the  bandwidth  over  the 
thrust  capacity  more  than  it  should.  The  use  of  a  ratio  was  an  arbitrary  choice,  other 
functions  could  well  be  the  solution  to  the  problem.  Table  4  summarizes  the  dual  bum 
maneuver  results. 
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Figure  16  Single  Bum  Strategy  Orbital  Radius  (2  km  Bandwidth,  40  N  Thrust) 


Figure  1 7  Single  Burn  Strategy  Orbital  Radius  (2  km  Bandwidth,  80  N  Thrust) 
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Figure  18  Single  Burn  Strategy  Orbital  Radius  (2  km  Bandwidth,  160  N  Thrust) 
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Figure  1 9  Single  Bum  Strategy  Orbital  Radius  (2  km  Bandwidth,  320  N  Thrust) 


Figure  20  Single  Burn  Strategy  Orbital  Radius  (2  km  Bandwidth,  640  N  Thrust) 


Figure  21  Single  Burn  Strategy  Orbital  Radius  (2  km  Bandwidth,  1280  N  Thrust) 


Figure  24  Single  Burn  Strategy  Orbital  Radius  (10  km  Bandwidth,  160  N  Thrust) 


Figure  25  Single  Burn  Strategy  Orbital  Radius  (10  km  Bandwidth,  320  N  Thrust) 
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Figure  36  Single  Burn  Strategy  Orbital  Radius  (75  km  Bandwidth,  160  N  Thrust) 
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Figure  37  Single  Burn  Strategy  Orbital  Radius  (75  km  Bandwidth,  320  N  Thrust) 
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Figure  41.  Apogee,  Radius,  and  Perigee  (25  km  Bandwidth,  40  N  Thrust) 
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Figure  43 .  Dual  Burn  Strategy  Orbital  Radius  (2  km  Bandwidth,  80  N  Thrust) 
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Figure  44  Dual  Burn  Strategy  Orbital  Radius  (2  km  Bandwidth,  160  N  Thrust) 


Figure  47.  Dual  Burn  Strategy  Orbital  Radius  (2  km  Bandwidth,  1280  N  Thrust) 
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Figure  53  Dual  Bum  Strategy  Orbital  Radius  (10  km  Bandwidth,  1280  N  Thrust) 
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Figure  62  Dual  Burn  Strategy  Orbital  Radius  ( 


Figure  64  Dual  Burn  Strategy  Orbital  Radius  (75  km  Bandwidth,  640  N  Thrust) 


Figure  65  Dual  Burn  Strategy  Orbital  Radius  (75  km  Bandwidth,  1280  N  Thrust) 
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2.  Orbital  Path 

The  orbital  trajectory  for  the  dual  bum  maneuver  was  significantly  smoother  than 
that  of  the  single  bum.  Raising  the  perigee  radius  with  the  additional  thmsting  sufficiently 
reduced  the  eccentricity  to  enable  a  near-circular  decay  pattern  to  evolve  in  most  cases  In 
those  where  the  bum  started  early  the  eccentricity  was  decreased  enough  to  prevent  the 
wide  range  of  apogee  and  perigee  experienced  in  the  single  bum  maneuver 


Bandwidth 

[km] 

Thmst  [N] 

40 

80 

160 

320 

640 

1,280 

2 

SAT 

SAT 

SAT 

SAT 

+10% 

UNSAT 

SAT 

SAT 

SAT 

SAT 

+2% 

+12% 

25 

+30% 

+24% 

SAT 

SAT 

SAT 

+4% 

75 

+10% 

SAT  ^ 

+3% 

+5% 

+3% 

+3% 

SAT  =  Satisfactorily  Maintains  Band 
UNSAT  =  Unsatsifactory 

+_%  =  _  %  Above  Maximum  Radius  Limit 
(%  of  Bandwidth) 

Table  4  Dual  Bum  Strategy  Band  Maintenance  Summary 


C.  PROPELLANT  CONSUMPTION 

Propellant  consumption  is  calculated  for  each  variation  of  the  single  and  dual  bum 
maneuvers,  as  well  as  for  an  FKT  at  the  maximum  radius,  midband  radius,  and  minimum 
radius  An  additional  consideration  involves  the  strategy  of  using  a  FKT  to  keep  the 
satellite  at  the  upper  radial  limit  until  the  time  when  removing  any  thrusting  would  cause 
the  orbit  to  decay  to  the  minimum  radius  by  the  end  of  the  mn  period.  TFCL  consumption 
is  compared  against  these  four  base  usage  patterns. 

Consumption  for  both  the  single  and  dual  burn  strategies  is  approximately  the 
same.  The  exact  value  at  any  given  instant  may  be  higher  for  one  or  the  other,  but  the 
average  values  overall  are  generally  equal.  Plateaus  between  burns  are  longer  in  the  dual 
bum  case,  as  expected.  Bums  occur  more  frequently  in  the  single  burn  maneuver  because 
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the  perigee  radius  may  not  be  sufficiently  raised  by  any  given  burn  to  preclude  drag 
quickly  decaying  it  below  the  minimum  limit. 

Comparison  of  the  two  methods  with  the  four  baselines  described  above  provides 
insight  into  their  efficiency.  In  all  cases  considered  the  propellant  required  by  either  TFCL 
strategy  is  slightly  less  than  or  roughly  equal  to  the  midband  FKT  counterpart.  As  the 
bandwidth  gets  very  large  the  TFCL  values  appear  to  drop  away  from  the  midband  ones. 
The  low  altitude  FKT  consumes  more,  significantly  more  in  cases  of  larger  bandwidths, 
than  the  midband  FKT  The  upper  radius  FKT  and  its  modified  version  use  less  that  the 
midband  FKT.  Both  TFCL  models  seem  to  be  bounded  above  by  the  midband  FKT  and 
below  by  the  high  altitude  FKT,  as  indicated  in  Figures  66  through  89.  This  suggests  that 
although  TFCL  is  better  than  attempting  an  FKT  at  the  middle  of  the  radial  bandwidth,  it 
is  less  efficient  than  either  an  FKT  at  the  maximum  radius  or  the  modified  maneuver 
described  previously. 
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Figure  66.  Propellant  Consumption  (2  km  Bandwidth,  40  N  Thrust) 
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Figure  67.  Propellant  Consumption  (2  km  Bandwidth,  80  N  Thrust) 


Figure  69  Propellant  Consumption  (2  km  Bandwidth,  320  N  Thrust) 
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Figure  70.  Propellant  Consumption  (2  km  Bandwidth,  640  N  Tluiist) 
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Figure  71.  Propellant  Consumption  (2  km  Bandwidth,  1280  N  Thrust) 
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Figure  72.  Propellant  Consumption  (10  km  Bandwidth,  40  N  Tlirust) 


Figure  73  Propellant  Consumption  (10  km  Bandwidth,  80  N  Tlirust) 
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Figure  74.  Propellant  Consumption  (10  km  Bandwidth,  160NTIirust) 


Figure  75.  Propellant  Consumption  ( 1 0  km  Bandwidth,  320  N  Thrust) 
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Figure  76.  Propellant  Consumption  (10  km  Bandwidth,  640  N  Tlirust) 
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Figure  77  Propellant  Consumption  (10  km  Bandwidth,  1280  N  Thrust) 
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Figure  78  Propellant  Consumption  (25  km  Bandwidth,  40  N  Tlirust) 
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Figure  79  Propellant  Consumption  (25  km  Bandwidth,  80  N  Thrust) 
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Figure  80  Propellant  Consumption  (25  km  Bandwidth,  160  N  Thrust) 


Figure  8 1  Propellant  Consumption  (25  km  Bandwidth,  320  N  Thrust) 
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Figure  82.  Propellant  Consumption  (25  km  Bandwidth,  640  N  Thrust) 
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Figure  83  Propellant  Consumption  (25  km  Bandwidth,  1280  N  Thrust) 
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Figure  84  Propellant  Consumption  (75  km  Bandwidth,  40  N  Thrust) 
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Figure  85  Propellant  Consumption  (75  km  Bandwidth,  80  N  Thrust) 
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Figure  S6.  Propellant  Consumption  (75  km  Bandwidth,  160  N  Thrust) 
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Figure  87  Propellant  Consumption  (75  km  Bandwidth,  320  N  Thrust) 
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Figure  88  Propellant  Consumption  (75  km  Bandwidth,  640  N  Thrust) 


Figure  89  Propellant  Consumption  (75  km  Bandwidth,  1280  N  Thrust) 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  purpose  of  this  thesis  was  to  develop  a  new  thruster  firing  control  logic  for 
either  a  single  burn  or  dual  bum  maneuver  that  successfully  maintained  an  orbital  band 
while  achieving  propellant  efficiency  approaching  or  exceeding  that  of  forced  Keplerian 
motion  Nondimensionalization  of  the  equations  permitted  more  accurate  computation 
and  easier  visualization  in  variation  of  parameters.  Results  indicate  that  both  strategies  are 
roughly  as  efficient  as  a  midband  FKT,  but  less  efficient  than  an  upper  limit  FKT  or  an 
upper  limit  FKT  with  decay  to  minimum  radius  at  end  of  life.  Reducing  the  eccentricity 
with  a  second  bum  increases  the  period  between  boosting  bums,  but  does  not  reduce 
propellant  consumption. 

The  question  of  optimality  of  thrust  cant  angle  is  opened  again,  based  on  the  new 
thruster  firing  control  logic.  Earlier  efforts  using  the  instantaneous  radius  and  osculating 
specific  energy  indicated  that  a  high  angle  was  required  to  keep  a  vehicle  within  a  radial 
band,  and  that  low  angles  were  incapable  of  constraining  the  radius.  Thmsting  along  the 
transverse  axis  using  the  osculating  perigee  and  apogee  radii  as  initiator  and  terminator, 
respectively,  is  capable  of  orbital  band  maintenance  and  significantly  more  efficient  than 
previous  methods 

In  a  separate  but  related  strategy  an  additional  bum,  based  on  the  flight  path  angle 
as  the  vehicle  approaches  apogee,  reduces  the  osculating  orbital  eccentricity,  allowing  a 
near-circular  decay  pattern  to  develop.  The  additional  thrusting  raises  the  perigee  radius, 
enabling  a  longer  time  between  orbit  boosting  maneuvers.  The  single  and  dual  bum 
maneuvers  are  equally  effective  at  band  maintenance  and  propellant  consumption,  but  the 
second  method  provides  a  potentially  more  desirable  orbital  pattern 

It  is  unclear  which  FKT  model  is  the  appropriate  baseline  to  compare  orbital  band 
maintenance  maneuvers  on.  At  the  low  altitude  of  the  band  the  FKT  solution  consumes 
more  propellant  because  it  is  fighting  a  denser  atmosphere.  Likewise,  the  high  altitude 
case  and  any  derivatives  from  it  have  the  benefit  of  having  a  significantly  rarefied  medium 
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to  travel  through,  and  require  less  propellant.  The  midband  version  is  not  the  most 
efficient,  thus  its  validity  as  the  basis  for  comparison  is  also  questionable 
The  results  of  this  thesis  indicate  several  areas  for  future  research: 

(1)  For  apogee  rasing  bums  fire  the  thruster  only  when  the  satellite  is  at 
perigee  The  initial  bum  in  both  methods  takes  place  when  the  osculating  perigee 
reaches  the  minimum  radial  limit,  regardless  of  the  position  of  the  vehicle  with 
respect  to  perigee.  Examining  at  each  occurring  perigee  the  radial  decay  from  the 
previous  perigee,  and  firing  the  thruster  if  the  drop  in  altitude  is  larger  than  the 
remaining  distance  to  the  lower  altitude  limit  could  reduce  propellant  consumption. 

(2)  Refining  the  thmster  initiation  for  the  second  bum.  Modification  of  the 
scaling  factor  or  exploration  of  other  controlling  functions  could  improve  the 
performance. 

(3)  Coupling  EITAG  with  a  control  logic  and  an  orbital  propagator.  The 
resulting  propellant  consumption  and  band  maintenance  could  then  be  compared 
with  FKT,  Pauls-Wilsey  and  this  method. 

(4)  The  effect  varying  the  specific  impulse  has  on  band  maintenance  and 
propellant  consumption. 

Finally,  additional  control  strategies  not  discussed  here  possibly  exist  that  could  provide 
the  answer  to  the  optimal  control  of  orbital  band  maintenance. 
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appendix 


A*  MAIN  PROGRAM  LISTING 

-  MAIN  PROGRAM 


Variable  Definitions 


a 

-  alphad 
^  alphar 
^  alto 

*  B 

-  Bbar 

*  beta 

*  D 

*  Dbar 

-  df 

*  dv 
^  dz 
i  DK 
%  e 

^  Esmax 
i  Esmin 
%  Es 
%  Eso 
^  Et 

^  g 

^  gamma d 

-  gamma r 
%  Gnot 

%  Gto 
%  h 
%  Isp 
%  L 

i  lambda 
%  Ifactor 
^  m 

*  mbar 
i  mf 

*  mfK 

-  mft 

^  mftK 
^  mnd 

*  mo 

-  orbits 

*  ra 

*  rho 

-  rhoalt 

*  rhofac 

*  rhonot 
^  rhostd 
t  Rmax 

-  Rmin 


=  Semi-Major  Axis 
=  Thrust  Angle  wrt  LH  [degrees] 

=  Thrust  Angle  wrt  LH  [radians] 

=  Midband  Altitude  [km] 

=  Ballistic  Coefficient  [kg/m^2] 

=  Nondimentionalized  Ballistic  Coefficient 
=  Atmospheric  Scale  Height 
=  Drag  [N] 

=  Nondimentionalized  Drag 
=  Drag  Factor 
=  Change  in  Velocity  [m/s] 

=  Orbital  Bandwidth  [km] 

=  FKM  Drag  [N] 

=  Eccentricity 

=  Maximum  Specific  Energy  at  Maximum  Radius  [J/kg] 

=  Minimum  Specific  Energy  at  Minimum  Radius  [J/kg] 

=  Specific  Energy  [J/kg] 

=  Initial  Specific  Energy  [J/kg] 

=  Total  Energy  [J] 

=  Standard  Earth  Gravity  [m/s^2] 

=  Flight  Path  Angle  wrt  LH  [degrees] 

=  Flight  Path  Angle  wrt  LH  [radians] 

=  Base  Flight  Path  Firing  Angle  [radians] 

=  Modified  Flight  Path  Firing  Angle  [radians] 

=  Angular  Momentum  [km'^2/s] 

=  Specific  Impulse  [s] 

=  Local  Computational  Variables 
=  Modified  Thrust  Firing  Scaling  Factor 
=  Length  Conversion  Factor  [1000  m/km] 

=  Spacecraft  Mass  [kg] 

=  Nondimentionalized  Mass 

=  Mass  of  Fuel  Burned  in  Time  Increment  [kg] 

=  FKM  Mass  of  Fuel  Burned  in  Time  Increment  [kg] 

=  Total  Mass  of  Fuel  Burned  [kg] 

=  FKM  Total  Mass  of  Fuel  Burned  [kg] 

=  Arbitrary  Mass  Used  for  Nondimentionalization  [kg] 
=  Initial  Spacecraft  Mass  [kg] 

=  Number  of  Orbits  Completed 
=  Apogee  Radius  [km] 

=  Calculated  Atmospheric  Density  [kg/m'^3] 

=  Reference  Atmospheric  Density  Altitude  [km] 

=  Density  Variation  Factor 
=  Reference  Atmospheric  Density  [kg/m^3] 

=  Standard  Earth  Density  [kg/m'^3] 

=  Maximum  Radius  of  Band  [km] 

=  Minimum  Radius  of  Band  [km] 
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r  Rminb 
r  ro 
i  rp 

*  t 

%  TBW 
^  Te 

*  tf 

-  tfl 

-  Th 

*  Thbar 

*  thetad 

*  thetar 
%  ThK 

-  Thin 

-  Thmbar 

*  tine 
^  tol 

^  Tpo 
%  Tp 

*  tstart 
%  tstop 

i  V 

*  Vbar 
Vimax 

*  Vimin 

*  Vmax 
^  Vmin 

r  Vprev 
%  x(l) 

%  x(2) 

%  X  (3) 

*  x(4) 

%  xbard) 

%  xbar{2) 

%  xbar(3) 

%  xbar{4} 

%  xbdot(l) 
%  xbdot{2) 
%  xbdot(3) 
r  xbdot{4) 


=  Nondimentionalized  Minimum  Radius  of  Band 
=  Initial  Orbital  Radius  [km] 

=  Perigee  Radius  [km] 

=  Time  [orbits] 

=  Base  Thrust-Bandwidth  Scale  Factor 
=  Earth's  Surface  Rotational  Period  [s] 

=  Final  Step  Time  [s] 

=  Thrust  Firing  Logic  Selector 
=  Thrust  [N] 

=  Nondimentionalized  Thrust 
=  Angle  from  Reference  Axis  [degrees] 

=  Angle  from  Reference  Axis  [radians] 

=  FKM  Thrust  [N] 

=  Maximum  (Blowdown)  Thrust  [N] 

=  Nondimentionalized  Maximum  (Blowdown)  Thrust 
=  Increment  of  Time  (Step  Size)  [s] 

=  Tollerance  Value  for  Computation 
=  Initial  Orbital  Period  [s] 

=  Orbital  Period  [s] 

=  Start  Time  [s] 

=  Stop  Time  [s] 

=  Velocity  [km/s] 

=  Nondimentionalized  Velocity 
=  Maximum  Velocity  at  Present  Orbit  [km/s] 

~  Minimum  Velocity  at  Present  Orbit  [km/s] 

=  Maximum  Velocity  at  Maximum  Radius  [km/s] 

=  Minimum  Velocity  at  Minimum  Radius  [km/s] 

=  Previous  Velocity  [km/s] 

=  Orbital  Radius  [km] 

=  Orbital  Radial  Velocity  [km/s] 

=  Orbital  Theta 

=  Orbital  Angular  Velocity  [1/s] 

=  Nondimentionalized  Orbital  Radius 

=  Nondimentionalized  Orbital  Radial  Velocity 

=  Nondimentionalized  Orbital  Theta 

=  Nondimentionalized  Orbital  Angular  Velocity 

=  Time  Derivative  of  xbar(l) 

=  Time  Derivative  of  xbar(2) 

=  Time  Derivative  of  xbar(3) 

=  Time  Derivative  of  xbar(4) 


clear  all 


^  Constants 


mu 

g 

Re 


=  3.98601208133E5; 
=  9.806; 

=  6378.145; 


%  Earth's  Gravitational  Parameter 
^  [km"‘3/s"2] 

%  Gravity  at  Earth's  Surface 
%  [m/s^2] 
i  Earth's  Radius 
i-  [km] 

%  Change  Degrees  to  Radians 
^  Convert  km  to  m 


dtor 

Ifactor 

tol 


=  pi/180; 
=  1000; 

=  le-6; 
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i count 


1:4; 


State  Variable  Numbers 


%  Independent  Variable  Initialization 


alto 

= 

260.0; 

dz 

= 

25.0; 

a 

= 

Re+alto+dz/2 ; 

alphad 

= 

0.0; 

B 

= 

150.0; 

betal 

46.9; 

df 

= 

1.0; 

dv 

o 

o 

e 

= 

o 

o 

o 

Isp 

300.0; 

index 

= 

1; 

mo 

= 

20000.0; 

mnd 

= 

20000.0; 

rhoalt 

= 

250.0; 

rhostd 

= 

7.248e-ll; 

rhof ac 

= 

12.47; 

rhonot 

= 

rhof a c* rhostd; 

Te 

= 

2*pi*  (  {  {Re'^3)  /mu)  "0.5)  ; 

thetad 

= 

180.0; 

Thm 

= 

320.0; 

Thf act 

= 

1.0; 

tine 

0.001; 

tstart 

o 

o 

tstop 

100; 

t 

tstart ; 

tfl 

= 

0; 

Gnot 

20e-6 

TBW 

=z 

10; 

lambda 

= 

TBW/  (Thm/dz''2)  ; 

Gto 

= 

lambda*Gnot; 

tol 

le-14; 

shoodi 

= 

0; 

doozit 

= 

100; 

%  Initializ 

:e  Orbital  Element  Variables 

thetar 

= 

thetad*dtor ; 

ra 

= 

a*  (14-e)  ; 

ro 

= 

{a* (l-e"2) ) / (l+e*cos (thetar) ) ; 

rp 

= 

a^*^  (1-e)  ; 

Tpo 

=r 

2*pi*sqrt ( (a"3) /mu) ; 

V 

= 

sqrt ( {2*mu/ro) - (mu/ a) ) ; 

ao 

= 

a; 

xd) 

= 

ro; 

x(2) 

= 

e^  (sin  (thetar)  )  *sqrt  (mu/  (a*  (l--e"2)  ) 

) ; 

x(3) 

thetar ; 

x(4)' 

= 

(sqrt  {a-*^mu'*^  (l-e"2)  ))/(x(l)"2); 

xbar ( 1 ) 

= 

X  ( 1 ) /ao; 
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xbar (2 ) 
xbar ( 3 ) 
xbar (4 ) 
Eso 
.  Et 
Es 


(Tpo^x (2 ) ) /ao; 
x(3)  ; 

(Tpo'^x(4)  )  ; 

(  {V"^2)  /2-inu/ro)  ; 
mo'*'Eso ; 


Initialize  Satellite  Variables 


m  =  mo ; 

alphar  =  alphad'^'dtor  / 

mbar  =  m/mnd; 

Thmbar  =  {Thm*Tpo^2) / (mnd*ao*lf actor ) ; 

Bbar  =  B/ ( rhonot*ao*lf actor ) ;  t  m/km 

Thbar  =  0.0; 

mft  =  0.0; 

mftK  =  0.0; 


^  Initialize  Altitude  Band  Constraint  Variables 

Rmax  =  ro+dz/2; 

Rmin  =  ro-dz/2; 

Rminb  =  Rmin/ao; 

Vmax  =  (mu/Rmax) ^0 . 5; 

Vmin  =  (mu/Rmin) ^0 . 5 ; 

Esmax  =  (  (Vmax'^2  ) /2  )  ~  (mu/Rmax)  ; 

Esmin  =  {  {Vmin'^2 ) /2 )  -  (mu/Rmin)  ; 

gamma  r  =  a  tan  (xbar  ( 2  )  /  (xbar  { 1 )  '^xbar  ( 4  )  )  )  ; 

norbs ( 1 )  =  0 ; 

Vimax(l)  =  Vmax; 


orad ( 1 ) 
oV(l) 
ogr ( 1 ) 
oh(l) 
oEs (1) 
oe(l) 
oa  ( 1) 
ora ( 1 ) 
orp (1) 
oTp(l) 
tf Im ( 1 ) 


V; 

0.0; 

ro"2'*^x(4}  ; 
Es ; 


Set  Forced-Keplerian  Motion  Drag 


rK{l) 

vK 


=  (ro-dz/2) /ao; 

=  Tpo*sqrt{ (2*rau/ (rK(l)*ao) ) - (mu/ (rK(l)*ao) ) )/a0; 


[rho,  Dbar]  =  leodrag(rK,  ao,  rhoalt.  Re,  beta,  vK,  Bbar,  df) 
DK  =  Dbar*mnd*ao*lf  actor/ {Tpo''2  )  ;  r  1000  m/km 
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^  Main  Program 


while  t  <=  tstop, 
while  index  <=  4, 

%  Determine  Nondimentionalized  Velocity 
[Vbar]  =  leovel {xbar ) ; 

*  Determine  Nondimentionalized  Drag 

[rho,  Dbar]  =  leodrag (xbar ,  ao,  rhoalt.  Re,  beta,  Vbar,  Bbar,  df) 

%  Determine  Local  Variables 

[L(l),  L(2),  L(3),  L{4),  L(5),  L(6),  L(7)]  =  leolocal ( xbar ,  mbar, 
Vbar,  Thbar,  Dbar,  alphar) ; 

%  Determine  Positional  Differential  Values 

[xbdot{l),  xbdot(2),  xbdot ( 3 ) ,  xbdot(4)]  =  leoxdot (xbar ,  L) ; 

%  Perform  RK4 

if  index  ==  1 

[xbt,  xbdt,  xbart]  =  leork4a (icount,  xbar,  xbdot,  tine) ; 
elseif  index  ==  2 

[xbdt,  xbart]  =  leork4b ( icount ,  xbt,  xbdt,  xbdot,  tine); 
elseif  index  —  3 

[xbdt,  xbart]  =  leork4c ( icount ,  xbt,  xbdt,  xbdot,  tine) ; 
else 

[xbart]  =  leork4d (icount,  xbt,  xbdt,  xbdot,  tine); 

end 

xbar (icount)  =  xbart ( icount ) ; 
if  rem( index,  2 )  ==  1 
t  =  ttO . S^tinc; 
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end 


index  =  index+1; 
end 

^  Reset  index  value 
index  =  1; 


^  Compute  Propellant  Consumption 

Th  =  Thbar*mnd'*'ao*lfactor/ (Tpo^2  )  ; 

mf  =  Th*tinc'*-Tpo/ (Isp*g)  ; 

m  =  m^mf ; 

mft  =  mft+mf; 

mfK  =  DK*tinc*Tpo/ {Isp*g) ; 

mftK  =  mftKtmfK; 

mbar  =  m/mnd; 


*  Increment  Output  Counter 
shoodi  =  shoodi+1; 

i  Store  Propellant  Consumption 


if  rem(shoodi,doozit)  ==  0 

mtv{fix (shoodi/doozit) +1)  =  mft; 

mtKv(fix (shoodi/doozit) tl)  =  mftK; 
mrat (fix (shoodi/doozit) +1)  =  mft/mftK; 
tflm ( fix ( shoodi/doozit) +1 ) =tf 1 ; 

end 


*  Determine  Orbit  Parameters 


eprev  =  e; 
grp=gammar ; 
Vprev  =  V; 

[r,  V,  gamma r^ 
[e,  a,  ra,  rp. 


h,  Es]  =  leoparml (xbar,  Tpo, 
Tp]  =  leoparm2(Es,  h^  mu); 


ao. 


mu) 


*  Thruster  Firing  Evaluation 
if  tfl  ==  2 

if  ra-r  <  r-rp 
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tfl  =  3; 


end 

end 

if  tlf  ==  3 

if  Thbar  >  0 

if  gammar-grp  >=  0 
tfl  =  1; 

end 

end 

end 

if  tfl  ==  0 

[Thval,  tfl]  =  leotflph(rp,  Rmin,  ra, 

elseif  rem  (tfl^2)  ==  1 

[Thval,  tfl]  =  leotflapn (gammar,  grp, 
Gto)  ; 

end 

Thbar  =  Thval'^Thf  act ; 

if  tfl  =-  4 
tfl  =  0; 

end 

if  e  <  tol 
e  =  0; 

end 

i  Output 

if  rein(shoodi,  doozit)  ==  0 

orad ( fix (shoodi/doozit ) +1)  =  r; 

oV(fix (shoodi/doozit) +1)  =  V; 

ogr ( fix ( shoodi/doozit ) +1 )  =  gammar; 

oh ( fix ( shoodi/doozit ) +1 )  =  h; 


Rmax,  Thbar,  Thmbar} ; 


V,  Thbar,  Thmbar,  tfl. 
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oEs ( f ix ( shoodi/doozit } +1 ) 

= 

Es  ; 

oe (fix (shoodi/doozit) tl) 

= 

e; 

oa ( fix ( shoodi/doozit ) +1 ) 

= 

a; 

ora ( fix ( shoodi/doozit ) +1 ) 

ra; 

orp (fix (shoodi/doozit) +1) 

= 

rp; 

oTp ( fix { shoodi/doozit) +1 ) 

= 

Tp; 

Thrst ( f ix ( shoodi/ doozit) +1 ) 

= 

Thbar ; 

oang ( f ix ( shoodi/doozit ) +1 ) 

*  Increment  Number  of  Orbits 

xbar ( 3 ) ; 

norbs ( f ix { shoodi/doozit } +1 )  =  (xbar ( 3 ) -x  ( 3 ) ) /  (2*pi )  ; 
end 

if  rem { shoodi , 10*doozit )  ==  0 
norbs { f ix ( shoodi/ doozit ) +1 ) 

end 

end 

end 

B.  VELOCITY  SUBROUTINE  LISTING 

-  NONDIMENTIONALIZED  VELOCITY  COMPUTATION 

*  Establish  as  Function 
function  [Vbarest]  =  leovel(xbar) 

^  Compute  Nondimentionalized  Velocity 

Vbarest  =  (  (xbar  (2  )  )  ^^2+  (xbar  (1)  *xbar  (4  )  )  ^2  )  *^0 . 5; 

return 

end 

C.  DRAG  SUBROUTINE  LISTING 

r  NONDIMENTIONALIZED  DRAG  COMPUTATION 
T  Establish  as  Function 

function  [rhogess,  Dbarest]  =  leodrag (xbar,  ao,  rhoalt.  Re,  beta,  Vbar, 
Bbar,  dfact) 

*  Determine  Density  Exponential  Term 
rhogess  =  exp (-( (xbar (1) *ao-Re) -rhoalt) /beta); 
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*  Compute  Nondimentionalized  Drag 

Dbarest  =  (df  act'*'rhogess*Vbar'"2  )  /  (2.0*Bbar)  ; 

return 

end 


D.  LOCAL  CONTRIBUTION  SUBROUTINE  LISTING 

*  LOCAL  CONTRIBUTION  SUBROUTINE 
%  Establish  as  Function 

function  [LI,  L2,  L3,  L4,  L5,  L6,  L7]  =  leolocal (xbar,  mbar,  Vbar, 
Thbar,  Dbar,  alphar) 

i  Define  Local  Computational  Variables 

LI  =  xbar (1) *xbar (4 ) *xbar (4) ; 

L2  =  4 . 0*pi*pi/ (xbar ( 1 ) *xbar ( 1) ) ; 

L3  =  (Dbar* (xbar (2) /Vbar) ) /mbar; 

L4  =  (Thbar*sin (alphar) ) /mbar; 

L5  =  (2.0*xbar(4)*xbar(2) )/xbar(l) ; 

L6  =  (Dbar/ (xbar (1) *mbar) )*( (xbar (1) *xbar (4) ) /Vbar) ; 

L7  =  (Thbar*cos (alphar) )/ (mbar*xbar (1) ) ; 

return 

end 


E.  DIFFERENTIAL  VALUES  SUBROUTINE  LISTING 

%  POSITIONAL  DELTA  VALUES  ROUTINE 
%  Establish  as  Function 

function  [dxl,  dx2,  dx3,  dx4]  =  leoxdot (xbar ,  L) 


*  Compute  Positional  Delta  Values 


dxl 

=  xbar (2)  ; 

dx2 

=  L(l)-L(2)-L(3)+L(4)  ; 

dx3 

=  xbar (4)  ; 

dx4 

=  -L(5)-L(6)+L(7)  ; 

return 

end 
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F. 


RUNGE-KUTTA  SUBROUTINE  LISTINGS 


i  RUNGE-KUTTA  4TH  ORDER  COMPUTATION 

*  Establish  as  Function 

function  [xbtemp,  xbdtemp,  nxbar]  =  leork4a (count,  txbar,  txbdot,  tminc) 

*  Perforin  first  part  of  RK4 

xbtemp (count)  =  txbar (count) ; 
xbdtemp ( count )  =  txbdot ( count ) ; 

nxbar ( count )  =  xbtemp ( count ) +0 . 5* tminc* txbdot ( count ) ; 

r  RUNGE-KUTTA  4TH  ORDER  COMPUTATION 

*  Establish  as  Function 

function  [xbdtemp,  nxbar]  =  leork4b (count,  txbar,  txbdot,  xbdot,  tminc) 

%  Perform  second  part  of  RK4 

xbdtemp (count)  =  txbdot (count) +2 . 0*xbdot (count) ; 
nxbar (count)  =  txbar (count) +0 . 5*tminc*xbdot (count) ; 

T  RUNGE-KUTTA  4TH  ORDER  COMPUTATION 
?  Establish  as  Function 

function  [xbdtemp,  nxbar]  =  leork4c (count,  txbar,  txbdot,  xbdot,  tminc) 

%  Perform  third  part  of  RK4 

xbdtemp ( count )  =  txbdot ( count ) +2 . 0*xbdot ( count ) ; 

nxbar(count)  =  txbar (count) +tminc*xbdot (count) ; 

*  RUNGE-KUTTA  4TH  ORDER  COMPUTATION 
i-  Establish  as  Function 

function  [nxbar]  =  leork4d ( count,  txbar,  txbdot,  xbdot,  tminc) 

*  Perform  fourth  part  of  RK4 

nxbar (count)  =  txbar (count) +tminc* (xbdot (count) +txbdot (count) ) /6; 

G.  ORBITAL  PARAMETER  SUBROUTINE  LISTINGS 

r  ORBITAL  PARAMETERS 
i  Define  as  Function 
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function [Pr,  PV,  Pgammar,  Ph,  PEs]  =  leoparml (xbar,  Tpo,  ao,  mu) 

*  Compute  Orbital  Parameters 
-  Radius 

Pr  =  xbar ( 1 ) *ao; 

*  Velocity 

PV  =  (ao*  (  (  (xbar  (2)  )  ^2+ (xbar  (1)  *xbar  (4)  ) ''2)  "0. 5)  ) /Tpo; 

^  Flight  Angle 

Pgammar  =  atan (xbar (2 ) / (xbar ( 1 ) *xbar { 4 ) ) ) ; 

r  Angular  Momentum 

Ph  =  Pr^PV^cos (Pgammar)  ; 

%  Specific  Energy 

PEs  =  (  (  (PV''2) /2)  -  (mu/Pr)  )  ; 

return 

end 

*  ORBITAL  PARAMETERS 
%  Define  as  Function 

function[Pe,  Pa,  Pra,  Prp,  PTp,  PVimin,  PVimax]  =  leoparm2(Es,  h 
%  Compute  Orbital  Parameters 

*  Eccentricity 

etest  =  1.0+2.0'^Es'*^  (  (h/mu)"^2)  ; 

Pe  =  sqrt ( etest ) ; 

-  Semi-major  Axis 

Pa  =  (-mu/ (2 , 0*Es)  )  ; 

*  Apogee  Radius 

Pra  =  Pa-*^  ( 1+Pe)  ; 

-  Perigee  Radius 

Prp  =  Pa^ (1-Pe) ; 


mu } 
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^  Period 


PTp  =  2 . 0*pi* (sqrt ( (Pa^3) /mu) ) ; 

*  Minmum  Velocity 

PVimin  =  {{mu^{l+Pe)}/(Pa^{l-Pe)))"0.5; 

*  Maximum  Velocity 

PVimax  =  (  (mu-*-  (l-Pe)  )  /  (Pa^  (1+Pe)  )  )  "0.5; 

return 

end 

H.  THRUSTER  CONTROL  SUBROUTINE  LISTINGS 

-  THRUSTER  FIRING  DETERMINATION  ROUTINE 

-  Establish  as  Function 

function  [Thrust,  tfl]  =  leotflph(Rp,  Rinin,  Ra,  Rmax,  Th,  Thmax) 
^  Thruster  Firing  Law 

*  Initialize  Thrust  Logic 
tl  =  0; 

^  Turn  Thrusters  Off 

if  Thcheck  ==  Thmax 

%  Check  Apogee  Radius 

if  Ra  >=  Rmax 

Thcheck  =  0.0; 
tl  =  2; 

end 

end 

-  Turn  Thrusters  On 
if  tl  ==  0 

if  Rp  <=  Rmin 

Thcheck  =  Thmax; 
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end 


end 

Thrust  =  Thcheck; 
tfl  =  tl; 
return 
end 

-  THRUSTER  FIRING  DETERMINATION  ROUTINE 
t  Establish  as  Function 

function  [Thrust,  tfl]  =  leotf lapn  (garnmar ,  grp,  V,  Thcheck,  Thmax,  tl, 
Gval) 

*  Thruster  Firing  Law 

*  Turn  Thrusters  Off 
if  tl  ==  1 

i  Check  Flight  Angle 

-  Turn  Thrusters  Off 

if  gamma r  >=  grp 

Thcheck  =0.0; 
tl  =  4;t0 

end 

end 

t  Turn  Thrusters  On 
if  tl  ==  3 

if  gamma r  >=0 

if  gamma r  <=  Gval 

Thcheck  =  Thmax; 

end 

end 

end 
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Thrust  =  Thcheck; 


tfl  =  tl; 


return 

end 
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